function [outputresultmap outputresulterrormap] = kneeregister(basedir,serie)
% kneeregister: run a dual rigid registration between different IR-SPGR
% inversion times
% pipeline: read in 2d dicom stack -> transform to 3d nifti -> draw square
% mask over bone -> determine lowest point of femur and highest point of
% tibia from cartilage masks -> perform separate registrations of femur and
% tibia on 3d Nifti images -> save back as 2d dicom stack

% Henk Smit

% Load image locations
cd(basedir)
load('locations.mat')

% Fixed image dir
IRdir = locations(serie).loc(1,1).dir;
IRfile = locations(serie).loc(1,1).file;
IRim = double(dicomread(fullfile(IRdir,IRfile)));
cd(IRdir);

% Volunteer serie dir
cd ..
volunteerdir=pwd;
d=dir;

% Postprocessing directory
cd ..
postprocdir = fullfile(pwd,'Postprocessing');
if ~exist(postprocdir, 'dir')   
      [success, message, messageid] = mkdir(postprocdir);
      if ~success    
         error( 'ratheartregister:invalidoutputdir', message );     
      end 
end

%Write dicom images in 1 3D nifti image
for i=3:size(d,1)
    if size(char(regexp(d(i).name,'SPGR', 'match')),2)>=1
        directory = d(i).name;
        [tstimg, info_structure] = loaddata([volunteerdir filesep directory filesep], 1 , 1.5);
        tstimg = squeeze(tstimg);
        T = affineTransformFromDicominfo( info_structure );
        [nii] = make_niiT(uint16(tstimg), T);
        save_nii( nii, [volunteerdir filesep directory, filesep,'3D.nii'])
    end
end

% Fixed image
fixedim=fullfile(IRdir, '3D.nii');

% Moving images
j=0; %counter
for i=3:size(d,1)
    if (size(char(regexp(d(i).name,'SPGR', 'match')),2)>=1 && size(char(regexp(d(i).name,'2100', 'match')),2)<1) %SPGR but not 2100 (fixed)
        j=j+1;
        directory = d(i).name;
        movingims(j).name=fullfile(volunteerdir, d(i).name, '3D.nii');
        movingims(j).dir = fullfile(volunteerdir, d(i).name);
    end
end

% find elastix.exe
currentdir= fileparts( which('dGEmricIRregister.m') );
cd(currentdir)
cd ..
elastixdir = [pwd filesep 'elastix'];

% A square around the knees for masks
screensize = get(0,'Screensize');
figure;
imshow(IRim,[1.1*min(IRim(:)) 0.5*max(IRim(:))],'InitialMagnification',100)
set(gcf, 'Position', [0.1*screensize(3),0.1*screensize(4),0.7*screensize(3),0.7*screensize(4)]);
title(['Click, hold and draw a square around the knee bone. Avoid the high intensities at the skin'])
RegMask=imrect;
position=wait(RegMask);
close;


%Read in femur_wb and tibia mask
    cd(fullfile(postprocdir,'Mapping_Masks'));
    dmask=dir;
for j=3:size(dmask,1)
    m=0;
    x=dmask(j).name;
    for mask=3:size(dmask,1) 
        if size(char(regexp(dmask(mask).name,'fem_wb', 'match')),2)>=1 && size(char(regexp(dmask(mask).name,'IR', 'match')),2)>=1
            load(dmask(mask).name,'M')
            femurMask=M;
        elseif size(char(regexp(dmask(mask).name,'plat', 'match')),2)>=1 && size(char(regexp(dmask(mask).name,'IR', 'match')),2)>=1
            load(dmask(mask).name,'M')
            platMask=M;
        end
    end
end
    
% Find lowest femur and highest tibia voxel and 5 voxels to make sure mask includes cartilage
[platx,platy] = find(platMask~=0);
[femx,femy] = find(femurMask~=0);
upperPlat = min(platx)-5;
lowerFem = max(femx)+5;
    

% Create a femur and a tibia mask
femMask = zeros(size(tstimg));
femMask(position(2):lowerFem,position(1):position(1)+position(3),:)=1;
tibMask = zeros(size(tstimg));
tibMask(upperPlat:position(2)+position(4),position(1):position(1)+position(3),:)=1;

% Write mask to NIFTI
femregmaskdir = fullfile(postprocdir, 'femregmask.nii');
tibregmaskdir = fullfile(postprocdir, 'tibregmask.nii');
[femMasknii] = make_niiT(uint16(femMask), T);
[tibMasknii] = make_niiT(uint16(tibMask), T);
save_nii(femMasknii, femregmaskdir)
save_nii(tibMasknii, tibregmaskdir)

% Parameter file
par = fullfile(elastixdir, 'parameters_RigidLMI.txt');

for q=1:2 % Perform registration separately for femur and tibia
    for r=1:size(movingims,2) % Perform registration for all moving images
        display(['Performing registration ',num2str(r+size(movingims,2)*(q-1)),'/',num2str(2*size(movingims,2))])

        % Output directory
        if q==1
            outdir = fullfile(postprocdir,'Registrations','femur',movingims(r).name(end-9:end-7));
            if ~exist(outdir, 'dir')   
                  [success, message, messageid] = mkdir(outdir);
                  if ~success    
                     error( 'ratheartregister:invalidoutputdir', message );     
                  end 
            end

            fixedmask = femregmaskdir;
            movingmask = femregmaskdir;
            
        elseif q==2
            outdir = fullfile(postprocdir,'Registrations','tibia',movingims(r).name(end-9:end-7));
            if ~exist(outdir, 'dir')   
                  [success, message, messageid] = mkdir(outdir);
                  if ~success    
                     error( 'ratheartregister:invalidoutputdir', message );     
                  end 
            end
            
            fixedmask = tibregmaskdir;
            movingmask = tibregmaskdir;
        end

        movingim = movingims(r).name;

        cmdline1 = ['"' fullfile(elastixdir, 'elastix') ...
        '" -f "' fixedim '" -m "' movingim ...
        '" -fMask "' fixedmask '" -mMask "' movingmask ...
        '" -p "' par '" -out "' outdir '"'];

        % Call Elastix and perform registration
        [status1,result1] = system(cmdline1);

        % Save 3D dicom into 2D slices
        newinfo=dicominfo([movingims(r).dir filesep '0006.dcm']); %read info from arbitrarily one of the images in the serie
        objT = affineTransformFromDicominfo( newinfo );
        [img] = dicomread([outdir, filesep, 'result.0.dcm']);
        img = squeeze(img);
        write_dicom3D([outdir filesep], img, newinfo, objT) 
    end
end
